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We study the convergence properties of our implementation of the moving punctures approach at 
very high resolutions for an equal-mass, non-spinning, black-hole binary. We find convergence of 
the Hamiltonian constraint on the horizons and the L2 norm of the Hamiltonian constraint in the 
bulk for sixth and eighth-order finite difference implementations. The momentum constraint is more 
sensitive, and its L2 norm shows clear convergence for a system with consistent sixth-order finite 
differencing, while the momentum and BSSN constraints on the horizons show convergence for both 
sixth and eighth-order systems. We analyze the gravitational waveform error from the late-inspiral, 
merger, and ringdown. We find that using several lower-order techniques for increasing the speed of 
numerical relativity simulations actually lead to apparently non-convergent errors. Even when using 
standard high-accuracy techniques, rather than seeing clean convergence, where the waveform phase 
is a monotonic function of grid resolution, we find that the phase tends to oscillate with resolution, 
possibly due to stochastic errors induced by grid refinement boundaries. Our results seem to indicate 
that one can obtain gravitational waveform phases to within 0.05 rad. (and possibly as small as 
0.015 rad.), while the amplitude error can be reduced to 0.1%. We then compare with the waveforms 
obtained using the cZ4 formalism. We find that the cZ4 waveforms have larger truncation errors for 
a given resolution, but the Richardson extrapolation phase of the cZ4 and BSSN waveforms agree 
to within 0.01 rad., even during the ringdown. 
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I. INTRODUCTION 

Numerical relativity (NR) has progressed rapidly since 
the breakthroughs of 2005 QH3] that allowed for the long- 
term evolution of black- hole binaries (BHBs). Among 
NR's significant achievements are its contributions to- 
wards the modeling of astrophysical gravitational wave 
sources that will be relevant for the first direct detection 
and parameter estimation by gravitational wave obser- 
vatories [H [5]. NR has also made contributions to the 
modeling of astrophysical sources, notably, the model- 
ing of the recoil kick imparted to the remnant BH from 
a BHB merger due to unequal masses [SHE], the remark- 
able discovery of unexpectedly large recoil velocities from 
the merger of certain spinning BHBs [ M2"oT, |21)H2T)] , and 
the application of the numerical techniques to combined 
systems of BHs and neutron stars [301 - 138] . More mathe- 
matical aspects of relativity have also recently been inves- 
tigated, including the evolution of N-black holes [3"9"H4"T] , 
the exploration of the no-hair theorem [32 0S] , and cos- 
mic [44] and topological censorship [45] , as well as BHBs 
in dimensions higher than four |46, 47] . The current state 
of the art simulations can simulate BHBs with mass ra- 
tios as small as q = 1/100 [33 133 and highly spinning 
BHBs with intrinsic spins a = Sn/Mfj up to (at least) 
0.97 [5D] [3T]. Currently these runs are very costly and it 
is hard to foresee the possibility of completely covering 
the parameter space densely enough for match filtering 
the data coming from advanced laser interferometric de- 
tectors by the time they become operational. 



To reduce the computational costs, several low- 
accuracy approximation are sometimes used. Among 
them are the techniques introduced in Ref. [52] where 
the number of buffer zones at AMR boundaries is reduced 
by lowering the order of finite differencing by successive 
orders near the AMR boundaries, the use of simple in- 
terpolations of spectral initial data rather than using the 
complete spectral expansion [53] , and copying the ini- 
tial data to the two past time levels for use in prolonga- 
tion at the initial timestep. All of these approximations 
proved to be useful for numerical simulations, but each 
one also has the side-effect of introducing a (hopefully) 
small 0{h) error. In this paper, we examine the effects 
of these approximations by performing high-resolution 
simulations of equal-mass, non-spinning BHBs, a prob- 
lem generally considered well under control. We show 
that a non-convergent error, that cannot be detected by 
simple means, is present when these techniques are used 
together. However, even when low-accuracy approxima- 
tions are eliminated, an apparently stochastic error in 
the waveform phase is still present that prevents us from 
seeing clean convergence of the waveform even at very 
high resolutions. We estimate this stochastic phase er- 
ror is controllable to within NINJA and NRAR accuracy 
requirements, but does make it very difficult to get an 
unambiguous measurement of the waveform phase and 
phase error. 
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TABLE I: Initial data parameters. All runs used the same 
parameters. The punctures are locates at f = ±(a:,0, 0), 
with momenta ±(p r ,pt,0), puncture mass parameters m p , 
and zero spin. The ADM mass is 1.00000003. We use 
A^coio. x A^coio. x Ncdlo. collocation points for the TwoPunc- 
TURES spectral solver. 



x 4.250 


m p 


0.4887922277 


Pr 


-0.00114088 


pt 0.11081837 


iVcolo. 


40 


M H 


0.50580 



II. NUMERICAL RELATIVITY TECHNIQUES 

To compute the numerical initial data, we use the 
puncture approach [S3] along with the TwoPUNC- 
TURES [S3] thorn (see Table [I] for the initial data parame- 
ter). In this approach the 3- metric on the initial slice has 
the form r y a i, — (ipBL + w) 4 (5 a fc, where ipBL is the Brill- 
Lindquist conformal factor, 5 a b is the Euclidean metric, 
and u is (at least) C 2 on the punctures. The extrinsic cur- 
vature is given by Kij — (iJj + u)~ 2 Kij, where is a su- 
perposition of Bowen-York solution |55] for BHs with spin 
S (here zero) and momentum p. The Brill-Lindquist con- 
formal factor is given by i/jbl = 1 +J27=i rn i/(^\^'~' ^Dj 
where n is the total number of "punctures" , is a pa- 
rameter (not the horizon mass), and is the coordi- 
nate location of puncture i. We evolve these BHB data- 
sets using the LazEv !56 implementation of the mov- 
ing puncture approach [3] with the conformal function 
W = s/x = cxp(— 2(f)) suggested by Ref. [S7] (where \ 
is the evolution variable introduced in [2]). The moving 
puncture approach is based on the Baumgartc-Shapiro- 
Shibata-Nakamura (BSSN, BSSN-NOK) formalism [SJ- 
160] , where the gauge and evolution variables are adapted 
such that the system is finite at the punctures. For the 
runs presented here, we use centered, eighth-order finite 
differencing in space [39] and a fourth-order Runge Kutta 
time integrator. (Note that we do not upwind the advec- 
tion terms.) For points near (but not on) the compu- 
tational domain boundary, we reduce the order of finite 
differencing such that the stencils fit in the computa- 
tional domain (i.e. reduce from eighth to sixth to fourth 
to second). 

On the boundary points themselves, we use radiative 
boundary conditions for all variables, which on the x = 
x m j n and x = x max faces, have the form 

d t f\ x = -v(-d x f+^^), (1) 
\x r J 

where the wavespeed v and asymptotic value of the func- 
tion are parameters. We set the wavespeed to 1 
for 'jijjAij, and (3 l , while we set the wavespeed to \[2 
for a,K, and W. We calculate d x f using a one-sided, 
second-order stencil. We set /oo = for all variables 
except 7jj, a and W, where we set /oo = 1. 

In the method of lines approach, for every evolved func- 
tion /, there is a corresponding function /rhSi where 



dtf = /rhs- The evolution code fills in all the RHS grid- 
functions (including on the boundaries). No additional 
boundary conditions (except symmetry boundary condi- 
tions) are applied to the evolution variables themselves. 

After each substep of the RK4 integration, we enforce 
the algebraic constraints, 7 = det(7y-) = 1 and A = 
"f lJ Aij = by renormalizing 7^ and subtracting the trace 
from Aij, i.e. 

1 

7« -> ^Hj, 

and 

Aij — > A^ — -jijA. 

Our code uses the Cactus/EinsteinToolkit [6lTf63] 
infrastructure. We use the Carpet [M] mesh refinement 
driver to provide a "moving boxes" style of mesh refine- 
ment. In this approach refined grids of fixed size are ar- 
ranged about the coordinate centers of both holes. The 
Carpet code then moves these fine grids about the com- 
putational domain by following the trajectories of the two 
BHs. 

We use AHFinderDirect [65] to locate apparent 
horizons. We measure the magnitude of the horizon 
spin using the Isolated Horizons algorithm detailed in 
Ref. [66]. Note that once we have the horizon spin, we 
can calculate the horizon mass via the Christodoulou for- 
mula 

m H = yjKr + S 2 H tt*Tn*„), (2) 

where ra m = \J Aj (16n) and A is the surface area of the 
horizon (we neglect the effects of momentum). 

For the gravitational waveform, we calculate the 
Newman-Penrose ^4 Weyl scalar using fourth-order cen- 
tered finite differencing and then decompose ifj^ at fixed 
radii using spin-weighted spherical harmonics, ip r = 
Ylz c t,m—2Ytm- We use a fourth-order accurate integra- 
tion to calculate the {I, m) modes of ip4- We also cal- 
culate the Hamiltonian (H), momentum (C*), and BSSN 
constraint (Q l ) violations using fourth-order centered fi- 
nite differencing. We use fourth-order finite differencing 
to reduce the computational cost of calculating analysis 
quantities. 

For our gauge conditions, we use a modified 1+log 
lapse and a modified Gamma-driver shift condition [2, 
E3EH], and an initial lapse a(t = 0) = 2/(l + tp% L ). The 
lapse and shift are evolved with 

(d t -^di)a = -2aK, (3a) 
d t p a = (3/4)f a -?#\ (3b) 

where we use r\ = 2 for all simulations presented below. 

For our tests of the fast, but low-accuracy techniques, 
we performed a set of 11 simulations. In all cases, the 
physical boundaries of our computational domains were 



located at 400Af in all coordinate directions, although we 
did use 7r-symmetry and z-reflection symmetry to reduce 
the computational domain by a factor of 4. We chose con- 
figurations with a base outer resolution of ho = 400M/70 
with 9 levels of refinement. The finest resolution was 
M/44.8. We then used global refinements of this base 
configuration with outer resolutions of hi = 400M/84, 
h 2 = 400M/100, and h 3 = 400M/120. The ratio be- 
tween the timestep and spatial resolution on the coarsest 
grid was set to 0.06125, the next two levels had 1/2 this 
timestep, the following three had 1/4, and then each ad- 
ditional level had time refinement of a factor of 2 with 
respect to the next coarsest grid. The ratio k of the 
timestep and spatial resolution [Courant factor (CFL)] 
on the finest timelevel was k = dt/h — 0.49. In addition, 
we performed simulations with all CFL factors reduced 
by factors of 2, 2V2, and 4. For the low accuracy study 
we used the reduction of order technique proposed in |52j 
to reduce the number of buffer zones at the refinement 
boundaries to six. This is accomplished by reducing the 
order of finite differencing from eighth to sixth, fourth, 
third, and second on the fourth point, third, and sec- 
ond points from the refinement, respectively. On the 
boundary points themselves, a simple copy from the past 
timelevel is used. After the internal timestepping is com- 
plete, prolongation is then used to overwrite the data in a 
region six gridpoints deep on each refinement boundary. 
To initialize the past timelevels of the initial slice, we used 
a simple copy and we used the interpolation option of 
TwoPunctures to quickly interpolate the spectral so- 
lution on the grid. 

We performed 25 high-accuracy simulations. For these 
simulations, we used nearly identical spatial grid struc- 
ture, but had identical CFL factors on each gridlevel. In 
addition, we used resolution h% • • • /13 above and a higher 
resolution /14 = 400Af / 144 (we also performed a single 
run with resolution h 5 = 400M/174). We used CFL fac- 
tors of k = 1/4, k = 1/(4^2), and k = 1/8. A CFL 
factor of k = 1/2 was not stable (due to low-resolution 
in the outer grid [69] ). These runs also differed from 
the low-accuracy runs in that a full complement of 16 
buffer zones were also used and no reduction of order 
near AMR boundaries was necessary. In addition we 
used the init_3_timelevels algorithm to initialize the 
grid, which properly initializes the past two timelevels 
of the initial data slice, and we used the computation- 
ally expensive evaluation option for TwoPunctures. 
The final two parameters that determine the accuracy 
of our simulations are the dissipation order and prolon- 
gation order. We found that fifth-order Kreiss-Oliger 
dissipation was most effective at suppressing the high- 
frequency noise that originates as reflections from AMR 
boundaries (see Fig.[l]). We similarly used fifth-order pro- 
longation operators. As an additional check, we compare 
all runs with two other systems, an eighth-order system 
with seventh-order dissipation and prolongation (we use 
seventh rather than ninth to reduce the number of re- 
quired buffer points) and a fully consistent sixth-order 




FIG. 1: A plot of the 5R{i/) 4 } on the xy plane at t = 110M. 
Note the interference pattern that develops as part of the wave 
reflects off the AMR boundary. Also note that the reflection 
originate in the buffer region and not on the AMR boundary 
itself. The reflected waves will soon interfere producing still 
higher spatial frequencies. 

system with seventh-order dissipation and prolongation 
operators. 

In the current work, we are interested in computing 
the phase of the dominant (t = 2, m = 2) mode of ^4 at 
finite radius. While there can be significant phase errors 
associated with extrapolation to r — 00, these errors are 
controllable by using larger computational domains, mul- 
tipatch methods [70], pseudospectral methods [71] . and 
Cauchy-Characteristic extraction 72, 73 . 

In the figures and tables below, we denote the low ac- 
curacy simulations by "L855" and higher-accuracy sim- 
ulations by "H855" , "H877" , and "H677" . Here the first 
digit indicates the spatial finite-difference order, the sec- 
ond indicates the dissipation order, and the third indi- 
cates the spatial prolongation order. 



III. RESULTS 

Our studies were motivated by phase errors in 
intermediate-mass-ratio simulations 74J. In that study, 
we found that mass-conservation was an important crite- 
rion for an accurate evolution. Here we study the conser- 
vation of the horizon mass as a function of resolution and 
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TABLE II: Differences between the low-accuracy, and high- 
accuracy configurations. 'N Buffer' is the number of buffer 
zones at refinement boundaries, 'Initial Data' refers to the 
method used for computing the conformal factor from the 
spectral coefficients, 'AMR Initialization' refers to the order 
of approximation used to initialize the past time levels of the 
initial slice (used for prolongation only). Here, 'radius' refers 
to the "half-width" in all directions of each level and CFL 
refers to the ratio between dt and h on that level. All high- 
accuracy simulations (H855, H877, H677) use the same grid 
structures and initialization parameters. 





low 


high 


N Buffer 


6 


16 


Initial Data 


interpolation 


evaluation 


AMR Initialization 


1st order 


4th order 





low 


high 


level 


radius 


CFL 


radius 


CFL 





400 


k/8 


400 


K 


1 


208 


k/8 


208 


K 


2 


115 


k/8 


115 


K 


3 


60 


k/4 


60 


K 


4 


30 


k/4 


30 


K 


5 


12 


k/2 


12 


K 


6 


5 


K 


5 


K 


7 


1.5 


K 


1.5 


K 


8 


0.75 


K 


0.75 


K 



100 200 300 

t/M 



400 



500 



Similarly, we can describe horizon fluxes of a vector field 
C l with the flux integral 



FIG. 2: The horizon mass as a function of time for different 
resolutions and CFL factors for the L855 algorithm (top) and 
H855 algorithm (bottom). Note the poor mass conservation 
for the low accuracy simulations. Also note that the mass for 
the H855 simulations is unaffected by CFL (i.e. the curves cor- 
responding to different CFLs lie on top of each other), while 
the mass is strongly affected by CFL for the low accuracy 
simulations. 



CFL factor for both the L855 and H855 configurations. 



A. Horizon Mass 

The horizon mass for non-spinning BHs is given by the 
irreducible mass M irr = -^4^/16, where 



s, 



(4) 



the integral is performed over the surface of the horizon 
and cPs is the proper volume element on the horizon as- 
sociated with the induced metric on the horizon. We can 
also define the horizon average of a function H as 



(H) = 



Hd 2 s 
U 2 s ' 



(5) 



(O = 



§C i R,d 2 



(6) 



where Ri is the unit norm on the horizon. 

As shown in Fig. [2j the conservation of the horizon 
mass for the L855 simulations is strongly affected by the 
CFL factor, with the largest deviations being between 
the k = 0.49 and K = 0.49/4 CFL factors. On the other 
hand, for the H855 simulations, the conservation of mass 
is unaffected by CFL, but is rather affected only by an 
increase in spatial resolution. Also note that the mass 
profiles are much flatter in the H855 simulations. 

We compare the horizon masses obtained with the 
H855, H877 and H677 systems (at a fixed CFL of 0.25) 
in Fig. [3j The H855 shows the most variation with res- 
olution, but appears to be comparable to H877 at high 
resolution. H855 and H877 both show a small drop in 
mass, which is unphysical, while H677 shows an increase 
in mass. At intermediate resolutions, it appears that 
H877 or H677 are better than H855 in terms of mass 
conservation. At high resolutions, H855 and H877 ap- 
pear to be slightly more accurate than H677 since they 
show flatter profiles of the mass versus time. 
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FIG. 3: The horizon mass versus resolution for the H855, 
H877, and H677 systems (all with CFL=0.25). Here the hori- 
zon mass has been smoothed by using a running average. At 
high resolutions, H855 and H677 both show very flat masses. 
H855 and H877 both show a small drop in mass, which is 
unphysical, while H677 shows an increase in mass. 



B. Constraints 

The size of the constraint violations is a measure of 
how closely the simulations obey the Einstein equations, 
and, in addition, violations of the constraint can lead 
to mass loss. We examined the Hamiltonian constraint 
violation both in the bulk and the average of the Hamil- 
tonian constraint violation on each horizon. The idea of 
using the horizon average constraints is to see how much 
strange matter is falling into the BHs (and thus a source 
of mass fluctuation). 

As seen in Figs. [4] and [6| the Hamiltonian constraint 
violation in the bulk is significantly larger for the L855 
simulations than the higher accuracy H855 simulations, 
with the constraint violation appearing immediately at 
the start of the simulation and decreasing towards the 
values seen in the H855 simulations. The H855 simula- 
tions, on the other hand, show clear fourth-order con- 
vergence. The constraints on the horizons converge to 
fourth-order for the H855 simulations, with similar val- 
ues to the L855. We can therefore conclude that ab- 
sorption of constraint violation by the horizon is not the 
ultimate source of either the lack of conservation in the 
mass or the phase error in the waveform. Spikes in the L2 
norm of the Hamiltonian, which are an artifact of the of 
7r-symmetry (possibly in the way the grids around each 
BH need to be enlarged near the symmetry boundary, or 
even an analysis artifact), as shown in Fig. [5] have been 
removed from Fig. [4] In Fig. [5j we show the the phase 
difference between a run that uses 7r-symmetry and one 
that does not. The differences in the phase are under 
0.001 rad. over the entire waveform. 

In Figures [7] and [8j we compare the L2 norm and 
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FIG. 4: Convergence of the L2 norm of the Hamiltonian con- 
straint in the volume between the horizons and a coordinate 
sphere at r = 20Af. Fourth-order convergence, consistent 
with the constraint calculation algorithm (and also the time 
integrator) is apparent for the H855 simulations. The low 
accuracy simulation shows a larger error during the inspiral, 
which may be the source of the error in the simulation. Spikes 
associated with the BHs crossing the symmetry boundary at 
x = 0, which are an artifact associated with -n symmetry, have 
been removed (see Fig. p5J|. 



horizon-averaged values of the Hamiltonian constraint vi- 
olations for the H855, H877, and H677 systems (all with 
CFL 0.25). Unlike the horizon mass, here H855 shows 
better behavior than H877. That is H855 shows conver- 
gence of the L2 norm to zero, while H877 shows con- 
vergence to a finite value. In Fig. [9j we show the L2 
norm of the x component of the momentum constraint 
(other components are qualitatively similar). Here, only 
the H677 simulations show consistent convergence with 
resolution. 

To summarize, at medium resolutions, both H877 and 
H677 preserve the horizon mass better than H855. At 
high resolution, H855 appears to be a little better than 
H877 or H677 (slightly flatter profile). The horizon- 
averaged constraint violations appear to be smallest for 
H677, followed by H855. The L2 norm of the con- 
straint violation favors H677 prior to merger (and af- 
ter t — 150M), and H855 (high resolution) post merger. 
Prior to t = 150M, high resolution H855 is better. 

In Fig. [TO] we plot the horizon average of the momen- 
tum constraint {C 1 )h and BSSN constraint (G l )h, where 
Q l = V + dj^'i , which is identically zero if the underly- 
ing ADM equations are solved exactly. These two hori- 
zon quantities can be thought of as the flux of constraint 
violation entering the horizons. For the momentum con- 
straints, the three systems seem to have the same level 
of constraint violation. On the other hand, for the BSSN 
constraint, the H855 system shows larger constraint vio- 
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FIG. 6: Convergence of the horizon averaged value of the 
Hamiltonian constraint violation for the H855 system. The 
horizon average constraint violation of H855 is consistent with 
that of L855 (after accounting for differences in resolution). 
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FIG. 5: A plot showing that the spikes in the Li norm of 
the Hamiltonian are an artifact of 7r-symmetry when the BHs 
cross the symmetry boundary. The plot shows the Li norm 
for identical runs, with the exception that one does not use 
7T symmetry. Outside the spikes themselves, the two runs 
agree. The bottom plot shows the phase differences between 
(otherwise identical) runs with and without 7r-symmetry. 
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FIG. 7: The Li norm of the Hamiltonian constraint for the 
H855, H877, and H677 systems. H877 shows clear conver- 
gence to a finite value at all resolutions. For H855, the con- 
straint decreases with resolution but appears to converge to a 
smaller finite value. H677 shows the best behavior, but also 
converges to (a still smaller) finite value. Post merger, H855 
gives the smallest violation (but only at the highest resolu- 
tion). At lower resolution, H855 gives the largest violation. 
Note that the curves have not been rescaled. 
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FIG. 8: The horizon average value of the Hamiltonian con- 
straint for the H855, H877, and H677 systems. In all cases, 
the CFL was set to 0.25. The constraint violation have been 
rescaled assuming fourth-order convergence. In all cases, con- 
vergence to zero is poorest at the highest resolutions. H677 
shows the most convergent behavior, while H877 shows the 
poorest. 
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FIG. 10: The horizon momentum constraint flux {C 1 )h and 
BSSN constraint flux {G 1 )h- Here H677 results are displayed 
using dot-dashed curves, H877 using dashed curves, and H855 
using solid curves. The subplots are stacked by resolution, 
with hi on top, followed by hz, The bottom subplot shows 
results for h,4, as well as /15 (thick solid line) for H677. Note 
that both the momentum and BSSN the constraints have been 
rescaled by a factor if (hi/h^) 4 (i.e. assuming fourth-order 
convergence). 



lations at all resolutions than the other two. 



FIG. 9: The L2 norm of the a;-component of the Momen- 
tum constraint for the H855, H877, and H677 systems. Only 
the H677 system shows consistent convergence to zero with 
resolution. Note that the constraints have not been rescaled. 



Based on these results, one may expect that high- 
resolution H677 to be the most accurate. The horizon- 
averaged Hamiltonian and momentum constraint (but 
not BSSN constraint) seem to indicate that the next most 
accurate simulation, at high resolution, is H855. While 
it is surprising that H877 does not perform as well as 
H677, this may be due to effects of reduced dissipation 
combined with the eighth-order algorithm. 
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FIG. 11: The magnitude and phase of the (£ — 2,m — 2) 
mode of ip4 extracted at r — 50 M. The magnitude is shown 
on a log scale, while the phase is shown on a linear scale. 
The exponential decay of the quasinormal ringdown is evident 
after t > 585M. 



C. Phase errors 

For the L855 configurations we saw clean convergence 
of the phase of the (I = 2, m = 2) mode of tp4 (see 
Fig. fn]| at the smallest CFL factor (0.125). While eight- 



phase error is relatively large. In Fig. 12 we show the 
convergence of the phase for both the L855 and H855 
simulations. Fourth-order convergence is apparent in the 
H855 simulations, but the magnitude of the phase error 
is smaller by a factor of ~ 6 compared to the L855 error. 



In Figs. 13 and 
2, TO 
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we show the phase of the (£ = 
2) mode of ip4 as a function of CFL factor n. 
Note that this means that only the timestep is refined, 
not the spatial resolution. Once might expect that the er- 
ror would decrease, or at least not increase, as dt is made 
smaller. Generally, this is the case, however, we note that 
the error can increase as dt is made smaller due to the ef- 
fects of dissipation. For example, for a simple first-order 
discretization for the advection equation, at fixed spatial 
resolution, the error is non-monotonic as a function of 



CFL. In Fig. 15 we show the convergence of the phase as 
a function of the CFL factor (at fixed spatial resolution) . 
The L855 simulations show clear first-order convergence 
at all spatial resolutions, while the H855 simulations ac- 
tually show an increase in error with CFL at higher res- 



order convergence in the phase is apparent, the actual 



olutions. Combining the information from Figs. 15 and 
[13} we conclude that the L855 system show a clear trend 
where the waveform at infinite resolution is a function 
of CFL factor. This means that at fixed CFL, which is 
how numerical convergence studies are performed, there 
is a non-convergent, error in the phase. Fortunately, the 
phase for the H855 simulations appears to be indepen- 
dent of CFL. 



r 



D. 



Comparing waveform phases using different 
dissipation and prolongation operators 



In Fig 16 we show the phase of the (I = 2, m = 2) 
mode of -04 near merger for the H855, H877, and H677 
systems. A priori, we might expect that the H877 sys- 
tems is the most accurate, and would therefore show the 
smallest amount of scatter and spurious oscillations of 
the phase with resolution. In fact, this system, and H677, 
shows significantly larger scatter than H855. A possible 
explanation is that the H877 and H677 systems are the 
least dissipative, and therefore most sensitive to high- 
frequency grid noise errors (these high frequency signals 
originate as interference signals from the reflection of the 
initial data pulse at the AMR boundaries). 

The H877 system does not appear to be consistent with 
either H855 or H677 if we assume that the error in the 
phases is given by the differences between the two highest 
resolutions. However, we note that in all the cases, the 
phase is not a monotonic function of the resolution, i.e. 
there are oscillations of the phase. The differences be- 
tween the phases for H855 and H877 and H855 and H677 
are better measures of the error. In this case, we find 
that the phase error is closer to 0.05. We plot the differ- 
ence of the highest resolution H855 simulation with the 



highest resolution H877 and H677 simulations in Fig. [17] 
Note that, based on the differences between the H855 and 
H677 simulations, we predict a phase error, during ring- 
down, of 0.015 rad. (0.045 if we use H877). Interestingly, 
this lower number is very close to the difference between 
the H855 Z4 simulations (described below) and the H855 
simulations (see Fig. 20). 

We also noticed differences in the phase between the 
ET_2011_10 (and later) version of Carpet (see [52]) and 
the ET_2011_05 (and previous) versions. In Fig. 18 we 
use "hg" and "git" to denote the newer and older ver- 
sions, respectively. As seen in Fig. [18] the differences 
converge away rapidly with resolution. The source of the 
deviation appears to be a difference in the prolongation 
operators that affects non-smooth data. The exponen- 
tially growing error in the phase apparent during the 
late inspiral is one of the main reasons why obtaining 
highly- accurate phases is so difficult. 

To conclude this section, we mention that the phase 
error seen in the waveform is not a result of the extraction 
technique, but is rather associated with phase error in 
the orbital motion itself. This is demonstrated in Fig. 19 
where oscillations in phase (with resolution) are apparent 
in both the orbital phase and waveform phase at the same 
time. 
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FIG. 12: Convergence of the phase of the (1 = 2, m = 2) mode of 1J14 for the L855 (left) and H855 (right) simulations. In all 
cases the CFL factor was 0.125 on the finest grid. Higher-order convergence is apparent for the L855 system. However, the 
solution itself has a non-convergent error. 
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FIG. 13: Convergence of the phase with changes to the CFL factor. The top left panel shows the convergence of the phase for 
the L855 system. Slightly better than first-order convergence is seen. The top right panel shows the convergence rate, at the 
same spatial resolution, for the H855 simulations. Note the phase difference are a factor of 10-100 smaller, but not convergent. 
The bottom left shows the same convergence plot, but for the next higher spatial resolutions. Clean 4th-order convergence is 
seen at late times, at early times the differences are consistent with zero. Finally, the bottom right panel shows the same plot, 
but for the highest resolutions runs. Here again, convergence is not observed with CFL. Also the phase differences with CFL 
are larger than for the medium resolution. 



10 




629.2 629.3 629.4 629.5 629.6 629.7 629.8 629.9 
t/M 

FIG. 14: The phase of the (£ = 2,m= 2) mode of during 
the ringdown phase. Shown are phases for the H855 and L855 
at three different CFL factors, as well as the Richardson ex- 
trapolation phase based on the highest resolution runs with 
the smallest CFL. Note how the phase converges to differ- 
ent values depending on CFL for the L855 system, but shows 
consistency between CFL values for the H855 algorithm. 



E. Analysis using the cZ4 Formalism 

The waveforms do not appear to converge in a con- 
sistent manner at high resolutions. In order to try to 



elucidate the origin of this problem, we repeated several 
of the simulations using the Conformal Z4 algorithm de- 
veloped in Alic et al |75j . For the parameters, we 
used Ki — 0.1, K2 = 0, and K3 = (Alic et al suggest 
setting K3 = 1/2, but our simulations were based on an 
earlier version of their paper) . We used the LazEv cod- 
ing infrastructure to generate the Z4 evolution code. In 
addition, we modified the gauge conditions to use Eq. (§ 
above (with T l replaced by the cZ4 variable T l ), rather 
than the gauge condition used in |75) (for the damping 
parameters, we used the values given in |75)L For these 
Z4 runs, we used the H855 system, that is, eighth-order 
centered finite differencing, fifth-order dissipation, fifth- 
order spatial prolongation, a full complement of 16 buffer 
zones, and fourth-order accurate initialization. As is ev- 
ident in Fig. [20j for a given spatial resolution, the Z4 
algorithm exhibits a larger truncation error compared to 
BSSN. While the Z4 runs show better than 4th-order 
convergence, with clean convergence at all resolution at- 
tempted, the size of the error makes comparison with 
BSSN difficult. That is, with coarser resolutions, the 
BSSN algorithm also appears to be cleanly convergent. 
Nevertheless, we use a Richardson extrapolation of the 
waveform phase for Z4 and compare this extrapolated 
phase with the BSSN phases, as shown in Fig. 



20 Here 



the BSSN phases are all much closer to the extrapolated 
Z4 phase than any of the individual Z4 phases. It would 
thus appear that the Z4 and BSSN results are consistent. 



F. Convergence of the Amplitude 

We measure the amplitude variation with resolution, 
relative to the highest resolution, of the (t = 2, m = 2) 
mode of ip4 



camp _ 



(IMM-IMMVIMM) 



for the k = 0.25,0.177,0.125 H855 configurations. We 
find that the deviations are smallest for k = 0.25, larger 
for k = 0.177, and largest for n — 0.125. Interestingly, 
the k = 0.177 results are more noisy than either of the 
other two. As shown in Fig. |21[ clear convergence is seen 
for the n = 0.125 configurations, while an oscillation is 
seen in the n = 0.25 results, which is consistent with 
our previous findings for the phase errors. That is, while 
cleaner convergence is seen with k = 0.125, this appears 
to be a consequence of the larger error associate with 
smaller CFLs. Thus, for high accuracy runs, one should 
apparently not use n — 0.125. The overall relative am- 
plitude error appears to be controllable to within 0.1%. 

We conclude this section by noting that even the low- 
est accuracy L855 simulation is accurate enough for many 
applications where extreme phase accuracy is not needed 
(e.g. recoil studies). In Fig. 22 we show the orbital trajec- 
tory for the lowest resolution L855 and highest resolution 



H677 simulations. The two orbital trajectories overlap 
throughout the entire simulation, even at merger, where 
the phase error increases to 0.15 rad. 
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FIG. 15: Plot showing the phase for a short duration during 
ringdown for simulations using the H855 (top) and L855 (bot- 
tom) algorithms. The plots are arranged according to CFL 
factor re with k = 1/4 on the left, re = l/(4y / 2) in the middle, 
and re = 1/8 on the right. In each case, the x and y axis are 
the same for each panel for a given algorithm. For the H855 
runs, note that re = 1/4 curves are not in convergence order, 
while re = 1/(472) are in convergence order but the curve sep- 
arations are not consistent with convergence, re = 1/8 shows 
convergence, but the extrapolation is about 0.13 rad from the 
high resolution run which is much larger than expected. The 
L855 show apparently better convergence properties but the 
phases are not convergent to each other. 
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FIG. 16: The phase of the [l = 2, m = 2) mode of ip/t near 
merger for the H855, H877, and H677 systems. Note how the 
more accurate H877 and H677 runs actually show significantly 
more scatter at lower resolutions than H855. For both H877 
and H855, the phases for the /13 and /14 resolution lie nearly 
on top of each other. Given the phase differences with the 
other resolutions, this is likely due to an oscillation (i.e., non- 
monotonic dependence) in the phase with resolution. 




FIG. 17: The difference between the phase of the (£ = 2, m = 
2) mode of ■04 as calculated using the H855 and H877, as well 
as the H855 and H677 systems. The vertical line corresponds 
to the time when the waveform frequency is Mw = 0.2. At 
this time, the phase error appears to be as small as 0.005 
(0.012, if we use the differences between H855 and H877) rad. 
Given the results apparent in Fig. |16| we estimate this to be 
a reasonable error of the phase error. 
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FIG. 18: The phase of the (£ = 2, m = 2) mode of ip4 near 
merger for the H855 systems with a CFL of k = 0.25 using 
both the git and hg versions of the Carpet AMR driver. 
The phase deviation between the git and hg versions increases 
exponentially with time, but decreases rapidly with increasing 
resolution. For comparison, the phase difference between the 
medium and high resolutions runs is also shown. 




504.8 504.9 505 505.1 505.2 

t/M 



FIG. 19: A comparison of the phase differences in both the 
trajectory and waveform between the two highest resolutions 
for the H677 system. The plot ends at the point where the 
orbital separation is 0.05A/. The agreement between the two 
curves indicates that the error in the phase of the waveform 
is associated with the error in the phase of the trajectory. 
The bottom plot shows the orbital and waveform phases in a 
small time interval for the different resolutions. Note how the 
oscillation (with resolution) in the phase of the waveform are 
matched by similar oscillations in the orbital phase. 
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FIG. 20: (top left) The phase convergence for the Z4 evolutions and a comparison with the Richardson extrapolated phase 
from the Z4 evolutions with the corresponding BSSN phases, (top right) The waveform phase during the late ringdown. The 
BSSN phases all lie between the highest resolution Z4 phases and the extrapolated phases. However, there also seems to be a 
trend towards moving to more negative phases as the CFL is reduced. This is observable in both BSSN and Z4. (bottom left) 
The Z/2 norm of the Hamiltonian constraint for the Z4 runs with 2 different CFLs. Convergence with resolution and consistency 
between CFLs is apparent at early and late times. However Inconsistencies between CFLs and non-convergent behavior is 
apparent in the later inspiral and merger phases. The oscillations in H'HI^ during the inspiral is consistent with the findings 
of [75]. (bottom right) The L2 norm of the Hamiltonian constraint for the comparable BSSN runs. Here improvements in the 
CFL are apparent when increasing resolution and decreasing the CFL factor. (Note that the spikes apparent in Fig. [5] have 
been removed for clarity). 
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FIG. 21: The error in the amplitude 5 l amp of the (£ = 2,m = 2) mode of tp 4 (<5f mp = (\ip4(hi) - I^KMVI^KM))- Clean 
convergence is seen for the k = 0.125 runs, while a change in sign in £3 is seen in the k, — 0.25 runs. The magnitude of £3 is 
roughly a factor of 3 smaller for the k = 0.25 runs. 
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FIG. 22: A comparison of the orbital trajectory for the low- 
est resolution L877 (low-accuracy) simulation and the highest 
resolution H677 simulation. The phase error at merger is 0.15 
rad, which corresponds to a waveform phase error of 0.30 rad. 
The actual post-merger phase difference is about 0.20 rad. 



15 



IV. DISCUSSION 

After removing various low-accuracy approximations, 
we were able to reduce variation in the horizon mass of 
the two BHs in a BHB to the level of ftr 6 A/. Wc find 
the Hamiltonian constraint violations (both the Li norm 
in the bulk, and horizon-averaged values) are convergent. 
The momentum and BSSN constraints also converge to 
fourth-order on the horizons. Convergence of the mo- 
mentum constraint in the bulk is less clear, with only 
H677 showing clean convergence. 

We find a stochastic error in the phase of the waveform 
that prevents us from seeing convergence at high resolu- 
tion. We speculate that this stochastic error is associ- 
ated with the high-frequency, unresolved, wave from the 
initial-data burst of radiation that reflects off the AMR 
boundaries, producing a complex, highly-variable inter- 
ference pattern on the grid (see Fig. [T]) . Without a clean 
convergence of the high-resolution runs, it is difficult to 
predict the actual error in the phase of the waveform at 
merger and ringdown. We note that the oscillations in 
the phase (with resolution) appear to be controllable to 
about the level of 0.015 rad., which is accurate enough for 
data analysis purposes (i.e. exceeds NRAR and NINJA 
accuracy goals). Whether the phase error really is this 
small depends on the source of the oscillations. If these 
are due to essentially stochastic errors introduced by grid 
reflections, then it seems reasonable to assume the phase 
error is controllable to 0.015 rad. However, if the source 
of the oscillation is actually lower-order errors becom- 
ing more dominant at higher resolutions, then we need 
far higher-resolution runs. Since the lowest-order error is 
the second-order prolongation (in time) error, one may 
need to perform a convergence study (at very high res- 
olutions) where the grids are refined by factors of two 
between resolutions (with coarsest resolution as fine, or 
finer than, our finest resolution). Such a study would be 
computationally very expensive, but may ultimately be 
necessary. 

One reason why grid refinement boundaries may be a 
significant issue has to do with reflection of the (rela- 
tively) high-frequency initial burst of radiation. By con- 
struction, the AMR grids are well adapted to evolving 
the relatively low-frequency waveform signal (with small- 
est period of about 10M), but the initial burst has fre- 
quencies a factor of four higher. To properly evolve this 
pulse, may require a factor of four increase in resolution, 
which would require a factor of 256 in computational re- 
sources. Fortunately, this signal is not physical, and one 
only needs to confirm that the effects on the rest of the 
waveform of poorly resolving this pulse are within data 
analysis tolerances. However, use of data with smaller 
spurious high-frequency wave content |76l I77j and more 
correct low-frequency wave content [78H50] should help 
reduce AMR reflections and the corresponding stochas- 
tic phase error. We expect to recover full convergence 
of the moving puncture formalism in the unigrid limit 
performing tests similar to those in Rcf. 81J. 



In order to obtain an independent measure of the phase 
error in our simulations, we compare the phase of the 
waveform obtained using both the BSSN and cZ4 for- 
malisms. We find agreement at the level of 0.01 ra- 
dian, but only after a Richardson extrapolation. Be- 
cause both codes used very similar technologies and the 
same AMR implementation, it would be very useful to 
compare the results from these simulations with other 
codes, analogous to what was done with the Samurai 
project [82] • Comparisons with other AMR-based codes, 
such as BAM [55], and pseudospectral codes, such as 
SPEC [H3IH1, would be especially useful. 

Our simulations show that the L855 system is inap- 
propriate for use where high phase-accuracy is needed. 
However, for situations were lower accuracy is accept- 
able, e.g. for recoil studies, it can be used. Our tests 
also show that, even when using high-accuracy methods 
with a state-of-the-art AMR code for numerical relativity, 
clean convergence remains elusive, but we have seen im- 
provements with the prolongation order in systems such 
as H677 and H877, with signs of consistency for H677. 
We note that [70] use the numerically more expensive 
H899 algorithm, which also includes upwinded finite dif- 
ferencing stencils. They found that upwinding increases 
the accuracy of the simulation. 



V. CONCLUSION 

We analyze the accuracy of a BHB simulation by exam- 
ining the preservation of the individual horizon masses, 
convergence of the constraints, and convergence of the 
magnitude and phase of the (£ = 2,m = 2) mode tpi 
extracted at r = 50M. As seen in Fig. [2] the fast, but 
low-accuracy approximations lead to poor conservation 
of the mass compared to the slower, but more accurate, 
techniques. For the high-accuracy techniques, we find 
convergence of the constraints to the expected order, as 
seen in Fig. [5] A residual, possibly stochastic, phase error 
is seen in the waveform itself. By comparing the wave- 
forms generated using different techniques and different 
evolution systems, as shown in Figs. [T7] and [20] we find 
consistency in the phase at the level of 0.05 rad. over the 
entire waveform. 
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